function [yLevel_cu,xLevel_cup,diagOut] = policyFunctionEPer(xLevel_cu,model,setupEPer)

%% Simulation using hybrid approach
ny              = setupEPer.ny;
nx              = setupEPer.nx;
mx              = setupEPer.mx;
nx1             = setupEPer.nx1;
Gxx             = 1/2*reshape(model.gxx,ny,nx^2);
Gxxx            = 1/6*reshape(model.gxxx,ny,nx^3);
G4x             = 1/24*reshape(model.g4x,ny,nx^4);
Gssxx           = 6/24*reshape(model.gssxx,ny,nx^2);        %scaling = nchoosek(4,2)/24
G5x             = 1/120*reshape(model.g5x,ny,nx^5);
Gssxxx          = 10/120*reshape(model.gssxxx,ny,nx^3);     %scaling = nchoosek(5,2)/120
Gsssxx          = 10/120*reshape(model.gsssxx,ny,nx^2);     %scaling = nchoosek(5,3)/120

Hxx             = 1/2*reshape(model.hxx,nx,nx^2);
Hxxx            = 1/6*reshape(model.hxxx,nx,nx^3);
H4x             = 1/24*reshape(model.h4x,nx,nx^4);
Hssxx           = 6/24*reshape(model.hssxx,nx,nx^2);        %scaling = nchoosek(4,2)/24
H5x             = 1/120*reshape(model.h5x,nx,nx^5);
Hssxxx          = 10/120*reshape(model.hssxxx,nx,nx^3);     %scaling = nchoosek(5,2)/120
Hsssxx          = 10/120*reshape(model.hsssxx,nx,nx^2);     %scaling = nchoosek(5,3)/120

%Given x_t: compute x1_t+1 and y_t (endogenous variables)
%Standard perturbation given certainty equivalence
tic
x_cu    = xLevel_cu-model.h0;
y_cu    = model.gx*x_cu;
x_cup   = model.hx*x_cu;
if model.orderApp > 1
    AA_cu = kron(x_cu,x_cu);
    y_cu  = y_cu  + Gxx*AA_cu;
    x_cup = x_cup + Hxx*AA_cu;
end
if model.orderApp > 2
    AAA_cu  = kron(AA_cu,x_cu);
    y_cu    = y_cu  + Gxxx*AAA_cu;
    x_cup   = x_cup + Hxxx*AAA_cu;
end
if model.orderApp > 3
    A4_cu   = kron(AAA_cu,x_cu);
    y_cu    = y_cu  + G4x*A4_cu;
    x_cup   = x_cup + H4x*A4_cu;
end
if model.orderApp > 4
    A5_cu   = kron(A4_cu,x_cu);
    y_cu    = y_cu  + G5x*A5_cu;
    x_cup   = x_cup + H5x*A5_cu;
end

% computing y_t+1=y_cup
y_cup   = model.gx*x_cup;
if model.orderApp > 1
    AA_cup   = kron(x_cup,x_cup);
    y_cup    = y_cup  + Gxx*AA_cup;
end
if model.orderApp > 2
    AAA_cup  = kron(AA_cup,x_cup);
    y_cup    = y_cup  + Gxxx*AAA_cup;
end
if model.orderApp > 3
    A4_cup   = kron(AAA_cup,x_cup);
    y_cup    = y_cup  + G4x*A4_cup;
end
if model.orderApp > 4
    A5_cup   = kron(A4_cup,x_cup);
    y_cup    = y_cup  + G5x*A5_cup;
end

% The residuals
ferror           = DSGEforesight_N([y_cu;x_cup(1:mx,1)],x_cu,y_cup,1,setupEPer);
residualSum      = max(abs(ferror));
if residualSum < setupEPer.residualMax
    %Standard perturbation
    x1_cup = x_cup(1:nx1,1);
    diagSim = [residualSum,residualSum,0,0];
else
    %Extended path
    [Z0,y_tN,N]           = getStartingValues(x_cu,setupEPer);
    [x1_cup,y_cu,diagSim] = solutionExtendedPath(Z0,x_cu,y_tN,N,setupEPer);
end
timeStep               = toc;
diagSim_Hybrid        = [diagSim timeStep];
%Adding derivatives with respect to sigma to correct for uncertainty
if model.orderApp > 1
    x_cup(1:nx1,1)   = x1_cup + 1/2*model.hss(1:nx1,1);
    y_cu             = y_cu   + 1/2*model.gss;
end
if model.orderApp > 2
    x_cup(1:nx1,1)   = x1_cup + 3/6*model.hssx(1:nx1,:)*x_cu + 1/6*model.hsss(1:nx1,:);
    y_cu             = y_cu   + 3/6*model.gssx*x_cu          + 1/6*model.gsss;
end
if model.orderApp > 3
    x_cup(1:nx1,1)   = x1_cup +  Hssxx(1:nx1,:)*AA_cu + 4/24*model.hsssx(1:nx1,:)*x_cu + 1/24*model.h4s(1:nx1,:);
    y_cu             = y_cu   +  Gssxx*AA_cu          + 4/24*model.gsssx*x_cu          + 1/24*model.g4s;
end
if model.orderApp > 4
    x_cup(1:nx1,1)   = x1_cup +  Hssxxx(1:nx1,:)*AAA_cu + Hsssxx(1:nx1,:)*AA_cu + 5/120*model.h4sx(1:nx1,:)*x_cu + 1/120*model.h5s(1:nx1,:);
    y_cu             = y_cu   +  Gssxxx*AAA_cu          + Gsssxx*AA_cu          + 5/120*model.g4sx*x_cu          + 1/120*model.g5s;
end

%Given x_t: computing x2_t+1 (Exogenous states)
x_cup(nx1+1:nx,1) = model.hx(nx1+1:nx,1:nx)*x_cu(1:nx,1);

% Adding the steady state level
yLevel_cu  = model.g0 + y_cu;
xLevel_cup = model.h0 + x_cup;

% Diagostic output
diagOut.maxAbsfVal   = diagSim_Hybrid(:,1);
diagOut.maxAbsfVal_1 = diagSim_Hybrid(:,2);
diagOut.iter         = diagSim_Hybrid(:,3);
diagOut.solver       = diagSim_Hybrid(:,4);
diagOut.timePerPeriod= diagSim_Hybrid(:,5);

end

